!********************************************************************************************************************
! Module for numerical derivatives.
!
! Declare:
!
! type(tDerivative) :: derivative
!
! call new(derivative, h, [params], [universe], [lps], [chi], [galaxyBin])
!
! You have to specify the initial step size, h. Typically, 0.2 to 0.3 times the evaluation value of the derivative
! variable is a good starting point. The speed and smoothness of the derivative are sensitive to this value, so
! I recommend timing and testing.
!
! [params], [universe], [lps], [chi], [galaxyBin] are optional arguments to allow real(D) parameters and other types
! to be passed into the function to be differentiated. In the function, they can be access via get methods, e.g.:
!
! lps = getLPS(derivative)
!
! Now evaluate the derivative of a function at x_0 by calling:
!
! call derivate(derivative, function, x_0, answer)
!
! The function must have the interface function(derivative, x) for x a real(D)
!
! You can do
!
! call checkConverged(derivative, errorMessage)
!
! See all the good examples in galaxyBin/galaxyBin.F90
!
!********************************************************************************************************************
module mDerivative

        use mPrecision
        use mObject

        implicit none

        interface new
                module procedure new_derivative
        end interface new

        interface getParam
                module procedure getParam_derivative
        end interface getParam

        interface getParams
                module procedure getParams_2
                module procedure getParams_3
                module procedure getParams_4
        end interface getParams

        interface getUniverse
                module procedure getUniverse_derivative
        end interface getUniverse

        interface getLPS
                module procedure getLPS_derivative
        end interface getLPS

        interface getChi
                module procedure getChi_derivative
        end interface getChi

        interface getGalaxyBin
                module procedure getGalaxyBin_derivative
        end interface getGalaxyBin

        interface converged
                module procedure converged_derivative
        end interface

        interface checkConverged
                module procedure checkConverged_derivative
        end interface

        private
        public :: new
        public :: getParam, getParams, getUniverse, getLPS, getChi, getGalaxyBin
        public :: error, converged, checkConverged
        public :: derivate

contains

        subroutine new_derivative(derivative, h, params, universe, lps, chi, galaxyBin)

                type(tDerivative), intent(out) :: derivative
                real(D), intent(in) :: h
                real(D), intent(in), optional :: params(:)
                type(tUniverse), intent(in), optional :: universe
                type(tLPS), intent(in), optional :: lps
                type(tChi), intent(in), optional :: chi
                type(tGalaxyBin), intent(in), optional :: galaxyBin

                derivative%h = h

                if( present(params) ) then
                        if( size(params) > DERIVATIVE_MAX_PARAMS ) then
!                                derivative%error = .true.
                                write(0, *) &
                                'ERROR: Too many parameters given to derivative'
                                stop
!                                return
                        else
                                derivative%params( 1:size(params) ) = params
                        end if
                end if

                if( present(universe) ) derivative%u = universe
                
                if( present(lps) ) derivative%lps = lps
                if( present(chi) ) derivative%chi = chi
                if( present(galaxyBin) ) derivative%galaxyBin = galaxyBin

        end subroutine new_derivative

        !******************************************************************
        ! Calculate derivative of func using Ridders method, which uses
        ! polynomial interpolation to extrapolate the derivative to h -> 0.
        ! See Numerical Recipes in Fortran 77 ch. 5.7.
        !******************************************************************
        subroutine derivate(derivative, func, x, ans)

                use mMisc

                type(tDerivative), intent(inout) :: derivative
                interface
                        function func(derivative, x)
                                use mPrecision
                                use mObject
                                real(D) :: func
                                type(tDerivative), intent(in) :: derivative
                                real(D), intent(in) :: x
                        end function func
                end interface

                real(D), intent(in) :: x
                real(D), intent(out) :: ans

                integer, parameter :: MAX_ITERS = 10
                real(D), parameter :: STEP_FACTOR = 1.4_D, &
                                      STEP_FACTOR2 = STEP_FACTOR ** 2, &
                                      ! If error after an iteration is ABORT_FACTOR worse than the best so far, then
                                      ! quit the calculation. This may mean that the initial h is too large.
                                      ABORT_FACTOR = 2.0_D

                real(D) :: h
                integer :: i, j, ierrmin

                ! Use polynomial interpolation based on Neville's algorithm
                ! (see Numerical Recipes Ch. 3.1).
                ! These will be the x values:
                real(D) :: xs(MAX_ITERS - 1), errs(MAX_ITERS - 1)
                ! This is Neville's tableau:
                real(D) :: tableau(MAX_ITERS, MAX_ITERS)

                h = derivative%h

                ! I will take the x values to be the *SQUARE* of the STEP_FACTOR
                ! That is, xs = (1.4)**2, (1.4)**4, etc., though we calculate
                ! derivatives with spacing equal to h/1.4, h/(1.4)**2, etc.
                ! I'm not sure why its good to square STEP_FACTOR, but I will
                ! trust the Numerical Recipes authors, for now.
                call geometricSeries(xs, STEP_FACTOR2, STEP_FACTOR2)

                tableau(1, 1) = ( func(derivative, x + h) &
                                - func(derivative, x - h) ) &
                                / (2.0_D * h)

                derivative%error = huge(1.0_D)

                do i = 2, MAX_ITERS
                        h = h / STEP_FACTOR
                        ! This is P_i of Numerical Recipes in Fortran 77
                        ! figure 3.1.2:
                        tableau(1, i) = ( func(derivative, x + h) &
                                        - func(derivative, x - h) ) &
                                        / (2.0_D * h)
!                        write(0, *) 'tableau:', tableau(1, i)
                        ! Perform extrapolations "rightward" in the
                        ! tableau:
                        tableau(2 : i, i) = ( tableau(1 : i - 1, i) &
                                              * xs(1 : i - 1) &
                                            - tableau(1 : i - 1, i - 1) ) &
                                            / ( xs(1 : i - 1) - 1.0_D )
                        ! Compare each new extrapolation to one order lower,
                        ! both at the present stepsize and the previous one,
                        ! keeping the larger of the two:
                        errs(1 : i - 1) = max(abs( tableau(2 : i, i) &
                                                 - tableau(1 : i - 1, i) ), &
                                              abs( tableau(2 : i, i) &
                                                 - tableau(1 : i - 1, i - 1) ))

                        ierrmin = iminloc( errs(1 : i - 1) )
                        if( errs(ierrmin) <= derivative%error ) then
                                derivative%error = errs(ierrmin)
                                ans = tableau(1 + ierrmin, i)
!                                write(0, *) i, h, ans
                        end if

                        if(abs( tableau(i, i) - tableau(i - 1, i - 1) ) >= ABORT_FACTOR * derivative%error) then
!                                write(0, *) 'WARNING: higher order much worse in derivate, iter:', i
                                ! Make err relative:
                                derivative%error = derivative%error / ans
!                                write(0, *) '1:', derivative%error
                                return
                        end if

                end do

                ! Make error relative:
                derivative%error = derivative%error / ans
!                write(0, *) '2:', derivative%error
 
       end subroutine derivate

        real(D) pure function error(derivative)

                type(tDerivative), intent(in) :: derivative

                error = derivative%error

        end function error

        logical pure function converged_derivative(derivative)

                type(tDerivative), intent(in) :: derivative

                real(D), parameter :: TOLERANCE = 0.04_D

                converged_derivative = (derivative%error <= TOLERANCE)

        end function converged_derivative

        subroutine checkConverged_derivative(derivative, errorMessage)

                type(tDerivative), intent(in) :: derivative
                character(len = *), intent(in) :: errorMessage

#ifndef IGNORE_DERIVATIVE_CONVERGENCE
                if( .not. converged(derivative) ) write(0, *) 'WARNING: ', errorMessage, '.'
#endif ! IGNORE_DERIVATIVE_CONVERGENCE

        end subroutine checkConverged_derivative

        real(D) pure function getParam_derivative(derivative)

                type(tDerivative), intent(in) :: derivative

                getParam_derivative = derivative%params(1)

        end function getParam_derivative

        pure subroutine getParams_2(derivative, p1, p2)

                type(tDerivative), intent(in) :: derivative
                real(D), intent(out) :: p1, p2

                p1 = derivative%params(1)
                p2 = derivative%params(2)

        end subroutine getParams_2

        pure subroutine getParams_3(derivative, p1, p2, p3)

                type(tDerivative), intent(in) :: derivative
                real(D), intent(out) :: p1, p2, p3

                p1 = derivative%params(1)
                p2 = derivative%params(2)
                p3 = derivative%params(3)

        end subroutine getParams_3

        subroutine getParams_4(derivative, p1, p2, p3, p4)

                type(tDerivative), intent(in) :: derivative
                real(D), intent(out) :: p1, p2, p3, p4

                p1 = derivative%params(1)
                p2 = derivative%params(2)
                p3 = derivative%params(3)
                p4 = derivative%params(4)

        end subroutine getParams_4

        type(tUniverse) function getUniverse_derivative(derivative)

                type(tDerivative), intent(in) :: derivative

                getUniverse_derivative = derivative%u

        end function getUniverse_derivative

        type(tLPS) function getLPS_derivative(derivative)

                type(tDerivative), intent(in) :: derivative

                getLPS_derivative = derivative%lps

        end function getLPS_derivative

        type(tChi) pure function getChi_derivative(derivative)

                type(tDerivative), intent(in) :: derivative

                getChi_derivative = derivative%chi

        end function getChi_derivative

        type(tGalaxyBin) pure function getGalaxyBin_derivative(derivative)

                type(tDerivative), intent(in) :: derivative

                getGalaxyBin_derivative = derivative%galaxyBin

        end function getGalaxyBin_derivative

end module mDerivative
